function project2()
close all
clear all


tic
[A Ak]= prelsi(40,80);
run_som(Ak,20,20,30000,'A80k40_')
run_som(A,20,20,30000,'A80_')

[A Ak]= prelsi(100,200);
run_som(Ak,20,20,30000,'A200k100_')
run_som(A,20,20,30000,'A200_')
toc


function [A Ak] = prelsi(k,p)
load A1000by501.mat

A0 = A1000by501;

A = zeros(p,501);
for i = 1:p
    r = 500+ceil(500*rand());
    A(i,:) = A0(r,:);
end

[m,n] = size(A);
[U,S,V] = svd(A);
N = min(m,n);
sigma = zeros(1,N);
for i = 1:N
    sigma(i) = S(i,i);
end
figure
semilogx(sigma)
title('singuler values')
saveas(gcf,strcat('A',num2str(m),'sigma.jpg'))
Uk = U(:,1:k);
Vk = V(:,1:k);
Sk = S(1:k,1:k);
Ak = Uk*Sk*Vk';


imagesc(log10(A))
colorbar
axis off
saveas(gcf,strcat('A',num2str(m),'a.jpg'))
spy(A)
axis off
saveas(gcf,strcat('A',num2str(m),'b.jpg'))
imagesc(log10(abs(Ak)))
colorbar
axis off
saveas(gcf,strcat('A',num2str(m),'k',num2str(k),'a.jpg'))
spy(Ak)
axis off
saveas(gcf,strcat('A',num2str(m),'k',num2str(k),'b.jpg'))


function run_som(Ak,m,n,step,name)

X0 = Ak;
[m1 n1] = size(X0);
X = zeros(m1,n1);
for term = 1:m1
    X(term,:) = X0(term,:)-mean(X0(term,:));
    X(term,:) = X(term,:)/max(X(term,:));
end

% Y = yi+1;

dstep = step/20;
mu = 0.1;
tau = step/2;
sigma0 = max(m,n)-2;
% w = SOM(X,Y,step,dstep,mu,n,m,tau,sigma0,'a.avi');
w = SOM(X,step,dstep,mu,n,m,tau,sigma0,name);



function w = SOM(X,step,dstep,mu0,n,m,tau,sigma0,aviname)
[N P] = size(X);
% w = 0.5*(rand(n*m,N) - 0.5);
w = (rand(n*m,N) - 0.5);
x = zeros(1,N);

fig = figure;
% aviobj = avifile(aviname);
% record = zeros(2,step);
for t = 1:step
    p = 1 + round((P-1)*rand());
    x = X(:,p)';
    i_c = argmin(x,w);
    for i = 1:n*m
        [h sigma] = neighborhood(i,i_c,m,n,t,tau,sigma0);
        mu = mu0*exp(-1*t/tau);
        w(i,:) = w(i,:) + mu*h*(x-w(i,:));
    end
    if mod(t,dstep)==0
        somedgeplot(w,m,n)
        title(strcat('step = ',num2str(t),' mu = ',num2str(mu),' sigma = ',num2str(sigma)))
        hold on
        umatplot(X,w,m,n)
        F = getframe(fig);
        saveas(gcf,strcat('som/',aviname,num2str(t),'.jpg'))
    end
end

function umatplot(X,w,m,n)
U = zeros(m,n);

for ix =1:m
    for iy = 1:n
        if ix==1
            dwr = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix+1,iy,m,n),:)).^2));
            dwl = dwr;
        elseif ix ==m
            dwl = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix-1,iy,m,n),:)).^2));
            dwr = dwl;
        else
            dwl = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix-1,iy,m,n),:)).^2));
            dwr = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix+1,iy,m,n),:)).^2));
        end
        if iy ==1
            dwd = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix,iy+1,m,n),:)).^2));
            dwu = dwd;
        elseif iy ==n
            dwu = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix,iy-1,m,n),:)).^2));
            dwd = dwu;
        else
            dwu = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix,iy-1,m,n),:)).^2));
            dwd = sqrt(sum((w(iPHI(ix,iy,m,n),:)-w(iPHI(ix,iy+1,m,n),:)).^2));
        end
        U(iy,ix) = mean([dwu,dwd,dwl,dwr]);
        x = 2.5+(ix-1)*5 + [-2 2 2 -2];
        y = -2.5 -(iy-1)*5 + [-2 -2 2 2];
        disp(U(iy,ix))
        fill(x,y,1-U(iy,ix))
    end
end



function somedgeplot(w,m,n)
[K,N] = size(w);
for ix = 1:n
    for iy = 1:(m-1)
        i1 = iPHI(ix,iy,m,n);
        i2 = iPHI(ix,iy+1,m,n);
        dw = sqrt(sum((w(i1,:) - w(i2,:)).^2));
        xplotfence(1-dw,ix,iy)
        hold on
    end
end

for iy = 1:m
    for ix = 1:(n-1)
        i1 = iPHI(ix,iy,m,n);
        i2 = iPHI(ix+1,iy,m,n);
        dw = sqrt(sum((w(i1,:) - w(i2,:)).^2));
        yplotfence(1-dw,ix,iy)
        hold on
    end
end
colormap('gray')
colorbar


function xplotfence(dw,ix,iy)
x = 2.5+(ix-1)*5 + [-1.5 -2.5 -1.5 1.5 2.5 1.5];
y = -5*iy + [-1 0 1 1 0 -1];
fill(x,y,dw)


function yplotfence(dw,ix,iy)
y = -2.5-(iy-1)*5 + [1.5 2.5 1.5 -1.5 -2.5 -1.5];
x = 5*ix + [-1 0 1 1 0 -1];
fill(x,y,dw)


function i_c = argmin(x,w)
[K,N] = size(w);
comp = zeros(K,1);
for i = 1:K
    comp(i) = sum((x - w(i,:)).^2);
end
[m i_c] = min(comp);


% the neighborhood function
% sigma should decrease as time propagates
function [h sigma] = neighborhood(i,i_c,m,n,t,tau,sigma0)
[ix iy] = PHI(i,m,n);
[i_cx i_cy] = PHI(i_c,m,n);
d2 = (ix - i_cx)^2 + (iy - i_cy)^2;
sigma = ceil(sigma0*exp(-t/(tau/log10(sigma0))));
h = exp(-1*d2/(2*sigma^2));


% PHI is the topology mapping from the neuron index to a 2D lattice
function [I,J] = PHI(i,m,n)
J = mod(i-1,n) + 1;
I = ((i-1)-mod(i-1,n))/m + 1;


% iPHI is the inverse of PHI
function i = iPHI(I,J,m,n)
i = (I-1)*n + J;



